Load network and scent data
From Tables S1 and S2 in: Kantsa, A., R.A. Raguso, A.G. Dyer, J.M. Olesen, T. Tscheulin, and T. Petanidou. 2018. Disentangling the role of floral sensory stimuli in pollination networks. Nature Communications 9: 1041. doi:[10.1038/s41467-018-03448-w](https://doi.org/10.1038/s41467-018-03448-w)
knitr::opts_chunk$set(echo = TRUE, comment=" ", cache=TRUE, fig.width=14, fig.asp=.678)
library(openxlsx)
library(data.table)
library(RColorBrewer)
pal <- brewer.pal(5,"Set2")
library(heatmaply)## Loading required package: plotly## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layout## Loading required package: viridis## Loading required package: viridisLite##
## ======================
## Welcome to heatmaply version 0.14.1
##
## Type citation('heatmaply') for how to cite the package.
## Type ?heatmaply for the main documentation.
##
## The github page is: https://github.com/talgalili/heatmaply/
## Please submit your suggestions and bug-reports at: https://github.com/talgalili/heatmaply/issues
## Or contact: <tal.galili@gmail.com>
## ======================library(bipartite)## Loading required package: vegan## Loading required package: permute## Loading required package: lattice## This is vegan 2.5-1## Loading required package: sna## Loading required package: statnet.common##
## Attaching package: 'statnet.common'## The following object is masked from 'package:base':
##
## order## Loading required package: network## network: Classes for Relational Data
## Version 1.13.0.1 created on 2015-08-31.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
## Mark S. Handcock, University of California -- Los Angeles
## David R. Hunter, Penn State University
## Martina Morris, University of Washington
## Skye Bender-deMoll, University of Washington
## For citation information, type citation("network").
## Type help("network-package") to get started.## sna: Tools for Social Network Analysis
## Version 2.4 created on 2016-07-23.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
## For citation information, type citation("sna").
## Type help(package="sna") to get started.## This is bipartite 2.08
## For latest changes see versionlog in ?"bipartite-package".
## For citation see: citation("bipartite").
## Have a nice time plotting and analysing two-mode networks.##
## Attaching package: 'bipartite'## The following object is masked from 'package:vegan':
##
## nullmodeldata(mosquin1967)
library(vegan)
library(phytools)## Loading required package: ape##
## Attaching package: 'ape'## The following object is masked from 'package:sna':
##
## consensus## Loading required package: mapslibrary(taxize)
#load network
net <- read.delim("Kantsa_etal_Supp_Data1.csv", skip = 2)
rownames(net) <- net[,1]
net[,1] <- NULL
net <- as.matrix(t(net))
snet <- sortweb(net) #sort by row & column totals
#load floral scents
scent.file <- "Kantsa_2018_Supp_Data2.xlsx"
plants <- getSheetNames(scent.file)
scentlist <- lapply(2:length(plants), read.xlsx, xlsxFile=scent.file)
plants <- plants[-1]
scent <- rbindlist(scentlist, fill=T, idcol="plant")[,-4]
colnames(scent)[2:3] <- c("cmpd","mean")
scent$plant <- factor(plants[scent$plant])
scent$cmpd <- factor(tolower(scent$cmpd))
cmpds <- levels(scent$cmpd)
scent.net <- dcast(scent, plant ~ cmpd, value.var="mean", fun.aggregate=mean, fill=0) #long to wide
scent.net <- as.data.frame(scent.net)
rownames(scent.net) <- scent.net[,1]
scent.net[,c(1,380)] <- NULL #weird NA column has one entry
scent.net <- as.matrix(scent.net)Plant phylogeny
#plants.class <- classification(plants, db="gbif")
#save(plants.class, file="plantsclass.Rdata")
load("plantsclass.Rdata")
plants.tree <- class2tree(plants.class)
plot(plants.tree)
Plant-pollinator network
Heatmap
heatmaply(net, scale="column")Matrix plot
#interaction matrix plot
visweb(net, type="nested", circles=T, boxes=F, circle.max=1.2)
Web plot
#interaction web plot
plotweb(net, text.rot=90, method="cca", col.high=pal[2], col.low=pal[1], col.interaction=pal[5])
Network topography
networklevel(net) connectance web asymmetry
0.06312657 0.63106796
links per species number of compartments
1.95631068 2.00000000
compartment diversity cluster coefficient
1.07901412 0.02631579
nestedness weighted nestedness
4.48059576 0.58672881
weighted NODF interaction strength asymmetry
10.80348739 0.20372363
specialisation asymmetry linkage density
-0.04424002 3.57255189
weighted connectance Fisher alpha
0.01734248 89.61252715
Shannon diversity interaction evenness
3.27629195 0.37393976
Alatalo interaction evenness H2
0.25285628 0.57919894
number.of.species.HL number.of.species.LL
168.00000000 38.00000000
mean.number.of.shared.partners.HL mean.number.of.shared.partners.LL
0.34944397 1.18776671
cluster.coefficient.HL cluster.coefficient.LL
0.34062967 0.21733898
weighted.cluster.coefficient.HL weighted.cluster.coefficient.LL
0.72421977 0.31277076
niche.overlap.HL niche.overlap.LL
0.12366249 0.14889538
togetherness.HL togetherness.LL
0.03019814 0.02302573
C.score.HL C.score.LL
0.75923950 0.74100046
V.ratio.HL V.ratio.LL
3.23242952 17.02718828
discrepancy.HL discrepancy.LL
248.00000000 250.00000000
extinction.slope.HL extinction.slope.LL
5.03260774 1.60178636
robustness.HL robustness.LL
0.81123903 0.61539429
functional.complementarity.HL functional.complementarity.LL
8044.31456516 7651.23911634
partner.diversity.HL partner.diversity.LL
0.94103312 1.28483306
generality.HL vulnerability.LL
2.67434342 4.47076037Scent network
Heatmap
heatmaply(scent.net, scale="column")Web plot
#interaction web plot
plotweb(scent.net, text.rot=90, method="cca", col.high=pal[2], col.low=pal[1], col.interaction=pal[5])
NMDS
scent.mds <- metaMDS(sqrt(scent.net), autotransform = F) Run 0 stress 0.2177614
Run 1 stress 0.2155518
... New best solution
... Procrustes: rmse 0.08183959 max resid 0.3861441
Run 2 stress 0.2277
Run 3 stress 0.2198936
Run 4 stress 0.2217743
Run 5 stress 0.2139972
... New best solution
... Procrustes: rmse 0.1078209 max resid 0.3550252
Run 6 stress 0.2201085
Run 7 stress 0.2306285
Run 8 stress 0.2158631
Run 9 stress 0.2417021
Run 10 stress 0.2264037
Run 11 stress 0.2330208
Run 12 stress 0.2175837
Run 13 stress 0.2506198
Run 14 stress 0.2138077
... New best solution
... Procrustes: rmse 0.07751348 max resid 0.2765495
Run 15 stress 0.2130163
... New best solution
... Procrustes: rmse 0.08566332 max resid 0.3435179
Run 16 stress 0.2426195
Run 17 stress 0.2136397
Run 18 stress 0.2156501
Run 19 stress 0.2146277
Run 20 stress 0.2169694
*** No convergence -- monoMDS stopping criteria:
20: stress ratio > sratmaxscent.op <- ordiplot(scent.mds, type="n")
points(scent.op, what="species", pch=3, col=pal[3])
text(scent.op, what="sites", cex=0.8)
Phylochemospace
scent.mds.taxa <- scent.mds$points
rownames(scent.mds.taxa) <- plants.tree$phylo$tip.label
plants.tree$phylo$edge.length <- plants.tree$phylo$edge.length+0.00001
phylomorphospace(plants.tree$phylo, scent.mds.taxa, control=list(col.node=0), label="none", lwd=1)
text(scent.mds.taxa[,1], scent.mds.taxa[,2], rownames(scent.mds.taxa), cex=0.8, offset=0.5, xpd=T)
Phylo Bar chart
par(mfrow=c(1,2))
par(mar=c(0.8,0,0.8,0))
plot(plants.tree, cex=1)
par(mar=c(0,0,0,0))
barplot(t(scent.net), col=sample(rainbow(ncol(scent.net))), names.arg=rep("",nrow(scent.net)), horiz=TRUE) 
par(mfrow=c(1,2))
par(mar=c(0.8,0,0.8,0))
plot(plants.tree, cex=1)
par(mar=c(0,0,0,0))
barplot(t(decostand(scent.net, method="tot")), col=sample(rainbow(ncol(scent.net))), names.arg=rep("",nrow(scent.net)), horiz=TRUE) 